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Molecular systematic investigation of Philippine puddle frogs 
(Anura: Dicroglossidae: Occidozyga Kuhl and Van Hasselt, 
1822) reveals new candidate species and a novel pattern of 


species dyads 


Kin Onn Chan", Sabine Schoppe?, Edmund Leo B. Rico” * And Rafe M. Brown? 


Abstract 


Focusing on the phylogenetic relationships of puddle frog populations spanning the biogeographic interface between 
Sundaland (Borneo) and the Philippines, we demonstrate, for the first time, a widespread geographic pattern involving the 
existence of multiple divergent and co-distributed (sympatric) evolutionary lineages, most of which are not each other’s 
closest relatives, and all of which we interpret as probable distinct species. This pattern of co-occurrence in the form of 
pairs of ecologically distinct puddle frog forms (dyads), prevails throughout northern Borneo, Palawan, Tawi-Tawi, the 
Sulu Archipelago, and western Mindanao (Zamboanga). Previously obscured by outdated taxonomy and logistical, legal, 
and security obstacles to field-based natural history studies, this pattern has remained hidden from biogeographers and 
amphibian biologists by an uncontested proposal that Philippine Occidozyga laevis is a single, “widespread,” and “highly 
variable” species. In this paper we use an integrative synthesis of new genetic data, organismal phenotypic data, historical 
literature reports, and ecological observations to elucidate an interesting and potentially widespread pattern of puddle frog 
species coexistence at the Sundaland—Philippine biogeographic interface. Calling attention to this pattern opens promising 
possibilities for future research aimed at understanding the scope of this dyads pattern, and whether it extends to the more 
northern reaches of the Philippines. On either side of Huxley’s and Wallace’s lines, data suggest that the majority of 
puddle frog dyads at a given locality are not each other’s closest relatives (are more distantly related, or non- 
monophyletic) and, thus, assembled ecologically, likely coexisting now as a result of their ecological tendencies toward 
distinct microhabitats (warmer stagnant pools in open areas, versus cool, flowing streams enclosed in forest). If these pairs 
of species types are determined to be the geographic norm among the more isolated, central, and northern, Philippine 
faunas, an obvious question will be whether they have evolved into dual ecological forms, possibly in response to 


ecological opportunity and/or reduced competition. 
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Introduction 


Puddle frogs of the genus Occidozyga Kuhl and Van 
Hasselt, 1822 are currently characterized with 12 species that 
are distributed throughout south and southeast Asia, ranging 
from eastern India through Bangladesh and southern China (as 
far east as eastern Fujian), and southwards throughout 
Indochina, the Thai-Malay peninsula, Indonesia (as far east as 
Sulawesi), Borneo, and the Philippines (Frost, 2020). Some 
species (referred to herein as “pond forms”) occur in muddy 
puddles or stagnant water bodies that are associated with 
disturbed habitats such as rice fields, road-side ditches, 
irrigation canals, and rubber plantations (Brown et al., 2013b; 
Devan-Song and Brown, 2012; Inger, 1954; Sanguila et al., 
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2016; Siler et al., 2011; Taylor, 1920), while others occur in 
swampy areas within secondary forests, ephemeral pools in 
primary forest, or stagnant, non-turbulent side pools 
disconnected from rivers (Brown et al., 2012; Inger, 1956, 
1954; Iskandar et al., 2011; Sanguila et al., 2016). However, 
some species such as O. baluensis (Borneo; Inger et al. 2017), 
O. semipalmata and O. tompotika (Sulawesi; Iskandar et al., 
2011) occur in shallow creeks with continuously running water. 
In the Philippines, one such species, O. diminutiva (Taylor, 
1922) has been reported exclusively from these forest- 
associated riparian microhabitats (Alcala and Brown, 1998; 
Brown and Alcala, 1970; Inger, 1954; Taylor, 1922). Most 
Occidozyga are morphometrically similar yet exhibit substantial 
intra- and interspecific variation in color-pattern among certain 
wide-ranging species complexes. To compound matters further, 
numerous species have overlapping ranges, resulting in 
taxonomic disarray stemming from the inconsistent and 
interchangeable application of different species names by 
different authors (AmphibiaWeb, 2019; Frost, 2020; Iskandar et 
al., 2011). 

Presently, two species of Occidozyga are known to occur 
in the Philippines: O. /aevis (Günther, 1858) and O. diminutiva 
(Taylor, 1922). Occidozyga laevis is known throughout the 
archipelago (Brown and Alcala, 1970; Inger, 1954; Taylor, 
1920), but the species’ type locality has only been reported as 
“Philippinen (=The Philippines)” without more specific island, 
or locality information, and this species also has been reported 
from Borneo (Inger et al., 2017). However, due to 
morphological similarities, O. laevis has long been confused 
with O. sumatrana on Borneo and other parts of Sundaland 
(Berry 1975; Boulenger 1882; Chan-ard 2003; Chan et al. 
2010a, b; Grismer et al. 2006; Inger 1966; Malkmus et al. 2002; 
Manthey & Grossmann 1997) and even now, their distributions 
are not satisfactorily characterized (Nutphund, 2001; Ohler, 
2003). As such, their validity as distinct species have not been 
adequately confirmed and O. laevis continues to be uncritically 
included in taxonomic and geographic summaries of Borneo's 
amphibian fauna (AmphibiaWeb, 2019; Frost, 2020; IUCN, 
2019) without accompanying data. 

In contrast, O. diminutiva is only known from three 
localities within the Philippines: the Zamboanga Peninsula of 
Mindanao Island, Basilan Island, and Jolo Island of the Sulu 
Archipelago (Brown and Alcala, 1970; Frost, 2020; Inger, 1954; 
Taylor, 1922) and since its description in 1922, no studies have 
critically assessed its validity other than Inger (1954) who 
evaluated its generic placement and transferred the species to 
Occidozyga from Micrixalus. Since that study, no additional 
comparisons or consideration of the species' taxonomic status 
and geographical distribution has been undertaken (Diesmos et 
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al., 2015). 

This study provides an initial molecular systematic 
investigation of Occidozyga in the southern Philippines with 
emphasis on populations that have had the potential of 
interacting with those from landmasses on the edge of the Sunda 
Shelf (Borneo Island). Given the biogeographical significance of 
this region (Brown et al., 2013a; Brown and Guttman, 2002; 
Brown and Alcala, 1970; Esselstyn et al., 2010; Lohman et al., 
2011), namely its role in contributing to the conceptualization of 
the discipline and the implications of Huxley's modification of 
Wallace's Line, which traditionally extended the edge of the 
Sundaic biogeographical faunal region to include Palawan 
(Inger, 1954; Mayr, 1944), focal studies specifically targeting 
fine scale differentiation and terrestrial faunal distribution 
patterns across this faunal zone interface have been surprisingly 
few (Brown, 2016; Brown and Guttman, 2002; Esselstyn et al., 
2010; Siler et al., 2012). 

Occidozyga constitutes an ideal study group by virtue of its 
distribution that spans Sundaland and the Philippines. We focus 
on northern Bornean and southern Philippine Occidozyga to 
examine the diversity and taxonomic status of these ecologically 
-variable, seldom studied, and potentially taxonomically 
confused amphibians. We ask three central questions: First (1), 
given that Occidozyga taxonomy and species names have been 
indiscriminately applied by various authors working exclusively 
from one side of the northern Malaysian-southern Philippine 
political boundary (e.g. O. sumatrana and O. laevis), can a 
survey of genetic data resolve taxonomic ambiguities and 
elucidate species diversity and distribution across this faunal 
zone interface at the edge of the Sunda Shelf? Second (2), given 
the traditional view of Palawan Island as a faunal zone extension 
of Sundaland, are current taxonomic arrangements tenable or 
reflective of evolutionary relationships? In other words, could 
Palawan support populations more closely-related to Sundaland 
species such as O. baluensis or O. sumatrana, especially given 
the latter’s notorious morphological similarity, which may 
render it indistinguishable from O. laevis? Finally (3) is it 
possible that the sympatric occurrence of species dyads (e.g. O. 
laevis and O. diminutiva, reported from the same areas but 
different microhabitats of the Sulu Archipelago and the 
Zamboanga Peninsula) could be more widespread? 


Materials and methods 


Morphological data and ecological observations 

Given the scope of this study and our goals of providing an 
initial survey of populations immediately on either side of the 
Sundaland-Philippine biogeographic interface (Fig. 1), we 
examined the same general organismal attributes that distinguish 
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Figure 1. Distribution of samples used in this study. Circles denote 
nominal taxa, whereas triangles represent candidate species. 
Phylogenetically unrelated (Fig. 2) and sympatric Puddle Frog species 
dyads are identified at three sites across the Sundaland-Philippine 
biogeographic interface, indicated by differently colored symbols 
(Borneo, Tawi-Tawi, Zamboanga); whereas genetically divergent but 
phylogenetically related (sister lineages) constitute the proposed dyad, 
identified on Palawan. 


the unambiguously distinct (valid) and sympatric species O. 
diminutiva and O. laevis: body size of males and females, 
general dorsal color pattern, ventral coloration, and 
microhabitat. Given that these categories of variation reliably 
allowed the last century of amphibian taxonomists to recognize 
these species (Alcala and Brown, 1998; Inger, 1954; Taylor, 
1922), we assessed the same variables to determine whether 
they could differentiate taxa on Palawan (a single species, O. 
laevis, most closely-related to nominal O. laevis from the 
oceanic islands of the Philippines), the Sulu Archipelago 
(represented in our data set by Tawi-Tawi, the closest 
Philippine island to Borneo, and where O. laevis has been 
reported), and western Mindanao (where the two species, O. 
diminutiva and O. laevis, have been reported but unconfirmed 
with genetic data or confirmations of the phenotypes or 
ecological attributes discussed above) For some populations 
(both forms/localities on Palawan, O. diminutiva and O. laevis 
in Zamboanga) we (RMB only, for precision) were able to 
collect body size data (snout-vent lengths in mm, of males with 
nuptial pads indicating sexual maturity) for 10-15 specimens 
using digital calipers; for Bornean O. cf. laevis and O. baluensis 
we augmented literature reports (Inger 1966; Inger et al. 2017) 
with confirmation of ranges, using additional measurements by 
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RMB from additional Bornean specimens (Table 2), to confirm 
comparability and repeatability of methodology between Inger’s 
(1966) and RMB’s data collection. At one locality (Tawi-Tawi), 
size data are unavailable (voucher specimens corresponding to 
our tissues destroyed; no other modern collections exist) and so 
the biologist who did the field work and observed these animals 
in life (ELBR) simply characterized these animals as “large” or 
“small,” based on memory, (such characterizations are included 
in quotes to indicate their subjectivity at this time). 


Molecular sampling and phylogenetic estimation 

For this study, we sequenced 38 samples from the holdings 
of the University of Kansas Biodiversity Institute (KU), and the 
Field Museum of Natural History, Chicago (FMNH). These 
include two outgroup samples (Hoplobatrachus rugulosus and 
Limnonectes woodworthii; Pyron and Wiens, 2011) and 36 
ingroup samples including Occidozyga sumatrana from 
Peninsular Malaysia and Indonesia; O. diminutiva from 
Zamboanga (type locality) and the Sulu Archipelago, 
Philippines; and O. laevis from Borneo, Palawan, and numerous 
other populations throughout the Philippines (Fig. 1). DNA was 
extracted with the Promega Maxwell© RSC Instrument using 
the Maxwell© RSC Tissue DNA Kit. We used the primers 16Sc 
-L (S-"GTRGGCCTAAAAGCAGCCAC-3^, and 16Sd-H (5'- 
CTCCGGTCTGAACTCAGATGACGTAG- 3’) to amplify the 
16S gene rRNA mitochondrial gene (Evans et al., 2003). 
Amplification was done using the following PCR thermal 
profile: 95 °C for 4min, followed by 35 cycles of 95 °C for 30 s, 
52 °C for 30 s, 72 °C for 70 s, and a final extension phase at 72 ° 
C for 7min (McLeod, 2010). Amplified DNA products were 
subsequently visualized on 1.0% agarose gels and sequenced at 
Genewiz, Frederick, MD. Sequences were assembled, aligned 
(MUSCLE algorithm), and concatenated in Geneious Pro 5.3 
(Kearse et al, 2012) prior to phylogenetic estimation. To 
provide a more comprehensive and accurate perspective of 
evolutionary relationships, we also supplemented our data with 
20 additional 16S sequences of O. lima and O. martensii from 
GenBank. All molecular samples used in this study and their 
associated GenBank accession numbers are listed in Table 1. 

We inferred phylogenies using maximum likelihood (ML) 
and Bayesian inference. The ML analysis was performed with 
IQ-TREE v1.6 (Nguyen et al, 2015) using the best-fit 
substitution model that was inferred by  ModelFinder 
(Kalyaanamoorthy et al., 2017). Branch support was calculated 
based on 1,000 bootstrap replicates using the ultrafast 
bootstrapping method (Hoang et al, 2017). A Bayesian 
phylogeny was inferred using BEAST v2.6 (Bouckaert et al., 
2019). The substitution model was estimated using bModelTest 
(Bouckaert and Drummond, 2017) and a relaxed log-normal 
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Table 1. Specimen identification, voucher information, localities, and Genbank accession numbers for sequences used in this study. 

















Taxa ie cnn RE GB# Museum Catalog # Locality 

From GenBank 

O. lima Java Island AB530619 O. lima AB530619 Indonesia: Java 

O. lima Vietnam AF206497 O. lima AF206497 Vietnam 

O. lima Thailand KR827958 O. lima KR827958 Thailand 

O. lima Cambodia KR827959 O. lima KR827959 Cambodia 

O. lima Laos KR827960 O. lima KR827960 Laos 

O. lima Myanmar MG935924 O. lima MG935924 Myanmar 

O. lima Myanmar MG935926 O. lima MG935926 Myanmar 

O. lima Myanmar MG935928 O. lima MG935928 Myanmar 

O. lima Kuala Lumpur, Malaysia AB488903 O. martensii AB488903 Malaysia: Kuala Lumpur 

O. martensii AF285214 O. martensii AF285214 Unknown 

O. martensii Vietnam DQ283357 O. martensii DQ283357 Vietnam 

O. martensii DQ458254 O. martensii DQ458254 Unknown 

O. martensii DQ458255 O. martensii DQ458255 Unknown 

O. martensii DQ458256 O. martensii DQ458256 Unknown 

O. martensii Ranong, Thailand AB530610 O. martensii AB530610 Thailand: Ranong 

O. martensii Thailand KP318725 O. martensii KP318725 Thailand 

O. martensii Myanmar MG935932 O. martensii MG935932 Myanmar 

O. martensii Myanmar MG935935 O. martensii MG935935 Myanmar 

O. martensii Myanmar MG935941 O. martensii MG935941 Myanmar 

O. baluensis Borneo, Malaysia DQ283143 O. baluensis DQ283143 FMNH 242747 Malaysia: Sabah, Borneo 

This study 

Hoplobatrachus rugulosus Hoplobatrachus rugulosus MT820164 UPLB(unknown) Philippines: Los Baños, 
Laguna, Luzon 

Limnonectes woodworthi Limnonectes woodworthi MT820165 KU 302234 Philippines: Catanduanes 

O. laevis Quezon, Luzon O. laevis MT820166 PNM (unknown) Philippines: Quezon, Luzon 

O. laevis Isabela, Luzon O. laevis MT820167 PNM (unknown) Philippines: Isabela, Luzon 

O. laevis Zamboanga KU 314470 O. laevis MT820168 KU 314470 Philippines: Pasonanca, 
Zamboanga, Mindanao 

O. laevis Zamboanga KU 314471 O. laevis MT820169 KU 314471 Philippines: Pasonanca, 
Zamboanga, Mindanao 

O. laevis Misamis Oriental KU 319796 O. laevis MT820170 KU 319796 Philippines: Misamis 
Oriental, Mindanao 

O. laevis South Cotabato O. laevis MT820171 PNM (unknown) Philippines: South Cotabato, 
Mindanao 

O. laevis Oriental Mindoro KU 302322 O. laevis MT820172 KU 302322 Philippines: Oriental 
Mindoro 

O. laevis Western Samar KU 306301 O. laevis MT820173 KU 306301 Philippines: Western Samar 

O. cf. diminutiva Tawi-Tawi ELR 161 O. cf. diminutiva MT820174 No voucher Philippines: Tawi-tawi 
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Table 1. Specimen identification, voucher information, localities, and Genbank accession numbers for sequences used in this study. 





Taxonomic identity 

















Taxa inferred in this study GB# Museum Catalog # Locality 

O. cf. diminutiva Tawi-Tawi ELR 170 O. cf. diminutiva MT820175 No voucher Philippines: Tawi-tawi 

O. cf. diminutiva Tawi-Tawi ELR 202 O. cf. diminutiva MT820176 No voucher Philippines: Tawi-tawi 

O. laevis Tawi-Tawi ELR 204 O. cf. laevis MT820177 No voucher Philippines: Tawi-tawi 

O. laevis Agusan Del Norte O. laevis MT820178 No voucher Philippines: Agusan Del 
Norte 

O. laevis Borneo FMNH 230140 O. cf. laevis MT820179 FMNH 230140 Malaysia: Sabah, Borneo 

O. laevis Borneo FMNH 230732 O. cf. laevis MT820180 FMNH 230732 Malaysia: Sabah, Borneo 

O. sumatrana Peninsular Malaysia FRIM 1132 O. sumatrana MT820181 FRIM 1132 Malaysia: Selangor 

O. sumatrana Peninsular Malaysia FRIM 1133 O. sumatrana MT820182 FRIM 1133 Malaysia: Selangor 

O. sumatrana Peninsular Malaysia FRIM 1936 O. sumatrana MT820183 FRIM 1936 Malaysia: Pahang 

O. sumatrana Java Island RMB 2132 O. sumatrana MT820184 MZB (unknown) Indonesia: Java Island 

O. sumatrana Java Island RMB 2133 O. sumatrana MT820185 MZB (unknown) Indonesia: Java Island 

O. sumatrana Java Island RMB 2134 O. sumatrana MT820186 MZB (unknown) Indonesia: Java Island 

O. cf. laevis 1 Irawan KU 326482 O. cf. laevis 1 MT820187 KU 326482 Philippines: Irawan, Palawan 

O. cf. laevis 1 Irawan KU 326483 O. cf. laevis 1 MT820188 KU 326483 Philippines: Irawan, Palawan 

O. cf. laevis 2 Brookes Point KU 326484 O. cf. laevis 2 MT820189 KU 326484 Philippines: Brookes Point, 
Palawan 

O. cf. laevis 2 Brookes Point KU 326485 O. cf. laevis 2 MT820190 KU 326485 Philippines: Brookes Point, 
Palawan 

O. cf. laevis 2 Brookes Point KU 326486 O. cf. laevis 2 MT820191 KU 326486 Philippines: Brookes Point, 
Palawan 

O. cf. laevis 2 Brookes Point KU 326487 O. cf. laevis 2 MT820192 KU 326487 Philippines: Brookes Point, 
Palawan 

O. cf. laevis 1 Irawan KU 308966 O. cf. laevis 1 MT820193 KU 308966 Philippines: Irawan, Palawan 

O. cf. laevis 2 Brookes Point KU 309476 O. cf. laevis 2 MT820194 KU 309476 Philippines: Brookes Point, 
Palawan 

O. cf. laevis 2 Brookes Point KU 309477 O. cf. laevis 2 MT820195 KU 309477 Philippines: Brookes Point, 
Palawan 

O. cf. laevis 2 Brookes Point KU 309478 O. cf. laevis 2 MT820196 KU 309478 Philippines: Brookes Point, 
Palawan 

O. cf. laevis 2 Brookes Point KU 309479 O. cf. laevis 2 MT820197 KU 309479 Philippines: Brookes Point, 
Palawan 

O. cf. laevis 1 Brookes Point KU 309480 O. cf. laevis 1 MT820198 KU 309480 Philippines: Brookes Point, 
Palawan 

O. diminutiva Zamboanga KU 321225 O. diminutiva MT820199 KU 321225 Philippines: Pasonanca, 
Zamboanga, Mindanao 

O. diminutiva Zamboanga KU 321226 O. diminutiva MT820200 KU 321226 Philippines: Pasonanca, 
Zamboanga, Mindanao 

O. diminutiva Zamboanga KU 321227 O. diminutiva MT820201 KU 321227 Philippines: Pasonanca, 





Zamboanga, Mindanao 
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model was used as the clock prior. Two independent MCMC 
chains were implemented (20 million generations each) and 
convergence was assessed using Tracer v 1.6 (Rambaut et al., 
2014). Converged trees were combined and a maximum clade 
credibility tree was inferred after the first 10% of sampled trees 
were discarded as burn-in. 


Candidate species delimitation 

We calculated uncorrected p-distances within and between 
nominal taxa using the complete-deletion function in MEGA-X 
(Kumar et al., 2018). To characterize and formally propose 
hypothesized (or “candidate”) species boundaries, we 
performed species delimitation analyses using mPTP (Kapli et 
al, 2017) and bGMYC (Reid and Carstens, 2012), both of 
which have been shown to be effective for single-locus datasets 
(Blair and Bryson, 2017; Dellicour and Flot, 2018; Tang et al., 
2014). For the mPTP analysis, we used the ML phylogeny as 
the input tree, and confidence of delimitation schemes were 
assessed using two independent MCMC chains at 5,000,000 
generations each. Support values represent the fraction of 
sampled delimitations in which a node was part of the 
Speciation process. To account for potential errors in 
phylogenetic estimation and uncertainty in model parameters, 
we implemented the bGMYC method that integrates over 
uncertainties in tree topology and branch lengths via MCMC 
(Reid and Carstens, 2012). As input, we used 100 randomly 
selected trees from the post burn-in combined runs of the 
BEAST analysis. For each tree, we ran the MCMC sampler for 
50,000 generations with a burn-in of 40,000, retaining 10,000 
post-burn-in generations with a thinning interval of 100. We 
then assessed species delimitation schemes from two 
conservative conspecificity probability thresholds (0.01 and 
0.05). 


Results 


Phylogenetic relationships 

Both ML and Bayesian phylogenies were concordant at all 
major nodes and all nominal species were monophyletic (Fig. 
2). Occidozyga lima and O. martensii were inferred as the first 
two branching lineages, respectively. Occidozyga baluensis was 
reciprocally monophyletic with O. diminutiva and this clade 
was the sister lineage to O. sumatrana (Fig. 2). Occidozyga 
laevis from the Philippine islands of Luzon, Mindoro, Samar, 
and northern Mindanao formed a highly structured clade, but 
populations from Palawan formed a clade that was reciprocally 
monophyletic with populations from Tawi-Tawi and Borneo. 
Within Palawan, two divergent and highly supported 
reciprocally monophyletic clades were inferred [Palawan form 1 
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(PA1) and Palawan form 2 (PA2); Fig. 2]. Branch support for O. 
laevis clades was high, while support for other species was 
mixed (Fig. 2). 


Candidate species delimitation 

Multiple highly divergent populations were detected within 
Occidozyga sumatrana, O. diminutiva, and O. laevis, indicating 
that distinct, undescribed species may potentially be present 
within each of these nominal species. Occidozyga sumatrana 
from Peninsular Malaysia and Indonesia (Java) were 8—9% 
divergent from each other, while O. diminutiva from Tawi-Tawi 
versus Zamboanga were 6-7% divergent (Fig. 3). Divergences 
within Occidozyga laevis were the widest, ranging up to 1596 
(Tawi-Tawi vs. Samar). Populations from Palawan were 6—9% 
divergent compared to populations from Borneo + Tawi-Tawi, 
whereas Palawan form 1 (river edges, puddles) and Palawan 
form 2 (small rapidly-cascading streams) were 5% divergent 
from each other. Occidozyga baluensis was 9-12% divergent 
compared to O. diminutiva. 

Both mPTP and bGMYC analyses supported the 
recognition of additional, undescribed species within 
Occidozyga diminutiva, O. sumatrana, and O. laevis. The mPTP 
analysis inferred O. diminutiva from Tawi-Tawi to be distinct 
from the Zamboanga population (Fig. 2). Occidozyga sumatrana 
was split into two species, represented by populations from 
Malaysia and Indonesia, while Occidozyga laevis was split into 
four species. At a conservative conspecificity probability 
threshold of 0.01, the bGMYC analysis lumped O. baluensis and 
O. diminutiva as a single species, but they were split at a 
threshold of 0.05. Similar to the mPTP analysis, Occidozyga 
sumatrana was also split into two species across both 
thresholds. At a threshold of 0.01, O. laevis was split into four 
species, while eight species were inferred at a threshold of 0.05 


(Fig. 2). 


Morphological attributes and microhabitat preferences 

Based on a combination of (1) existing taxonomy (Taylor 
1922; Inger 1954, 1966; Inger et al. 2017), (2) our assessment of 
traditional taxonomic character state differences from new 
localities and populations previously unreported in the literature 
(included herein), (3) available specimen-associated natural 
history museum data, field notes, new island records, and other 
occurrence records (summarized here), (4) our own field 
observations (microhabitat data, summarized together in this 
paper) and (5) limited body size information available from 
specimen data (Table 2) and the literature (Taylor 1922; Inger 
1954, 1966; Inger et al. 2017), we found a general pattern of two 
apparent forms (Table 2) of puddle frogs at each of the 
landmasses spanning the Sundaland-Philippine faunal zone 
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Table 2. Species identification for four proposed puddle frog species dyads, at four sites (Fig. 1), general subjective (putative) classification as 
“pond” versus “stream” ecological types (see text for discussion), microhabitats, forest cover, general size, coloration, and phylogenetic relation- 
ships (polyphyletic vs monophyletic; see Fig. 2) inferred for each dyad. Data sources and voucher specimen information are included for refer- 
ence. Table entries in quotes indicate qualitative impressions of field biologists in cases where no data are available and questions marks (?) indi- 
cate unknown. Data for Bornean species from Inger (1966) and Inger et al. (2017); Philippine populations data from Taylor (1920, 1922), Inger 
(1954) Alcala and Brown (1998), Brown et al. (2012, 2013b), Devan-Song and Brown (2012), Sanguila et al. (2016), Siler et al. (2011) with addi- 
tional size ranges augmented and these earlier reports confirmed with data collected from voucher specimens (by RMB), corresponding to the 
same individuals and populations from which genetic data presented here (Figs. 2, 3) were derived; all other size data from specimens deposited at 
KU (Table 2). Data sources: T22 (Taylor 1922); AB98 (Alcala and Brown 1998) I54 (Inger 1954); 166 (Inger 1966); I17 (Inger et al. 2017); T 
(This study). All specimens correspond to museum-deposited voucher specimens and genetic material (KU: https://collections.biodiversity.ku.edu/ 
KUHerps/; FMNH: https://collectionszoology.fieldmuseum.org/), unless one or the other has been lost/destroyed (voucher specimens from Tawi- 
Tawi Island [for which genetic material is available]). See Table 1 for correspondence between individual specimens from which genetic sequence 
data were derived. 
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s I Forest A Ventral Habitat Dyad 
Site/Taxa ecological female size à ; " Source Voucher 
type cover (SVL) coloration substrate relationship 
Zamboanga & Basilan 
O. cf. laevis pond open areas & 27-38; yellow muddy water polyphyletic T22, KU334623-29, 
forest 31-55 AB98, 334630-35, 
154,T 334636-37, 
33463841; 
O. diminutiva stream forest 14-24; Blotched stream bed T22, KU314464—69, 
22-29 white + I54, T 321225-53; 
brown AS 62559, 
CAS20069-71, 
0156-60. 
Tawi-Tawi & Jolo 
O. cf laevis pond open areas & "large" yellow muddy water polyphyletic T ELR204 (genetic 
forest sample) 
O. cf diminutiva pond + disturbed smaller" Blotched stream bed T ELR161, 170, 202; 
stream? forest white + CAS60681—83, 
brown 62528-34. 
Borneo 
O. cf laevis pond forest 21-31; yellow muddy water polyphyletic 166, FMNH29081, 
35—48 117, T 129091-111, 
129113, 129120; 
KU155621-22, 
15562324. 
O. baluensis stream forest 15-25; blotched stream bed I66 KU155617-20, 
25-35 white + 117, T 339545, 339547; 
brown FMNH36030-31, 
137444, 138075—84, 
138086-89, 
138096-101. 
Palawan 
O. cf laevis 1 pond open areas & 22-33; yellow muddy water monophyletic I54, T FMNH123585, 
forest 28-44 KU79025-32, 
91305, 309484—85, 
326482-83, 308966. 
O. cf laevis 2 stream forest 22-29.7; gray with stream bed T KU 79033, 
25.5-33.4 dark flecks 308966, 309164, 
and blotches 309476—80, 326482, 


326484-87. 
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Figure 2. Bayesian phylogeny based on 900 base pairs of the 16S mitochondrial rRNA gene juxtaposed with results from the mPTP and bGMYC 
species delimitation analyses. Gray bars represent species delimitation results that contradict current taxonomic arrangements. Support values at 
nodes are Bayesian posterior probabilities, followed by ML bootstrap support, and support values for the mPTP species delimitation analysis. 
Localities of interest are listed as BO = Borneo, TT = Tawi-Tawi, ZA = Zamboanga, PM = Peninsular Malaysia, IN = Indonesia, PA1 = Palawan 
form 1, and PA2 = Palawan form 2. Asterisks (*) denote locations where species dyads occur (BO, TT, and ZA). Inset photos of representative taxa 
include Occidozyga lima (I), O. martensii (II), O. diminutiva (III), O. sumatrana (IV), O. laevis nominal type from Philippines (V), and O. laevis 
from Borneo (VI). All photographs were taken by the authors except I, II, and IV, which were taken by Lee Grismer. 


interface  [Palawan, 


Tawi-Tawi, 


southwestern Mindanao 


(Zamboanga); each with two apparent morphological types], 
and Borneo (three or more forms). We referred loosely to these 
individual forms, which illustrate the duality of natural history 
observed at any single locality, as the “pond” form and “stream” 
form with the caveat that strong statements regarding adherence 
to this distinction and degree to which these forms actually 
exhibit truly more, or less specialization each will require 
verification (with substantially more data from detailed 
ecological studies and more sophisticated morphometric 
analyses based on large sample sizes of properly preserved 
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specimens), but based on the following general set of 
observations. Typically, at a given site or general area, with 
sufficient survey and re-survey work, field biologists have 
recorded two forms of puddle frogs: one form ("pond") 
constitutes a larger-bodied type, usually seen/observed/collected 
in more open and/or disturbed habitats, and which is 
encountered most often submerged in muddy puddles, stagnant 
small ponds, water buffalo wallows, flooded rice fields and 
irrigation conduits, or unconnected side-pools of slow-flowing, 
low-elevation, silted rivers near coastal areas. Another form (the 
“stream specialist”) has been detected by biologists in latter 
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stages of sampling, or has been detected only upon follow-up re 
-surveys in a particular area (Tawi-Tawi, Zamboanga) or 
following years of studies in an area (Palawan), perhaps 
suggesting differences in detectability which might reflect 
behavior and/or other ecological attributes and preferences—or 
just may reflect taxonomic or search-image bias on the part of 
field biologists. The “stream” form is usually smaller, it may be 
distinguishable with darker ventral coloration, also somewhat 
camouflaged dorsally, with spots and cross bars on the head, 
body, and limbs. In addition to its small body, darker, more 
cryptic coloration, the “stream” form appears associated with 
more enclosed, forested habitats, in some instances at higher 
elevation, and usually is not encountered in stagnant pools but, 
rather, in cooler forested areas, and situated in flowing or 
cascading streams. 

Available evidence suggests that the loosely-applied 
characterization (above) is not perfect and that some exceptions 
should be enumerated. Thus, on Languyan Island, Tawi-Tawi, a 
“pond” form was found in rice fields and buffalo wallows and a 
“stream” form (O. cf. diminutiva) was found in ephemeral 
puddles on the floor of a selectively-logged, secondary forest; 
other specimens were taken from a small stream (ELBR, pers. 
obs.). At Pasonanca, Zamboanga Peninsula of western 
Mindanao, a “pond” form was found primarily in muddy pools 
on the edge of an agricultural area, and submerged in water of 
unconnected side-pools of a large low-elevation river with, and 
a “stream” form (O. diminutiva) was encountered in flowing 
tributaries and small seeps, both flowing into the same river, but 
set back in the forest at steeper terrains where water was 
cascading over rocks and cliff faces (RMB, pers. obs.). On 
Palawan, the situation is less clear. In some southern areas in 
the vicinity of Brookes Point, Quezon, and Narra, only one 
microhabitat specialist, the “pond” form, has yet been recorded 
(both in stagnant pools in open areas and forest rivers and 
streams (RMB pers. obs.), whereas at others (Montible, Roxas, 
Taytay, San Vicente; SS and RMB, pers. obs.), a second form is 
suggested by our genetic data, smaller specimens in collections 
which were only collected in forested streams (RMB, pers. 
obs), and a tendency for specimens encountered in small 
forested streams to be more colorfully- and cryptically-blotched 
(SS, pers obs.). 


The puddle frog dyads pattern 

Connecting our disparate lines of evidence, even while 
acknowledging each’s associated caveats (e.g., the fact that our 
observations are from notes, records of field biologists, and 
specimen-associated data assembled by different investigators 
over a span of ~20 years), nevertheless led to the apparent 
conclusion of a dyads pattern of occurrence in puddle frogs 
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forms. The point at which our data intersected is the site-by-site 
(Fig. 1) accounting of our combined sources of information. 
Whether species dyads involved uncontroversial, named (Taylor 
1922; Inger 1954, 1966; Alcala and Brown 1998; Inger et al. 
2017), phylogenetically unrelated pairs of species (Borneo, 
southwestern Mindanao), a named species paired with a 
currently unrecognized candidate lineage (Tawi-Tawi), or a pair 
of haplotype clades which may correspond to disparate 
ecological observations and color pattern tendencies (Palawan: 
where the two forms may be sister species), an over-arching 
pattern of the apparently repeated duality of ecological and 
morphological types of puddle frogs exists at multiple localities 
immediately spanning Wallace’s and Huxley’s lines (Fig. 1). 
This is particularly evident when we considered that at most 
sites (1) our molecular phylogenetic analysis found strong 
Statistical support for the placement of two separate forms 
(unrelated in Borneo, western Mindanao, and Tawi-Tawi, but 
possibly sister lineages on Palawan); that in all known dyads, 
both forms are (2) genetically divergent lineages; and that (3) 
our DNA-facilitated identification of formerly recognized 
species (previously named species, now newly confirmed with 
genetic material from each named species’ type localities—and 
which enabled us to identify related but genetically distinct new 
allopatric candidate species such as O. cf. diminutiva on Tawi- 
Tawi, and O. cf. laevis in Borneo) suggests that other, putatively 
new, possibly undescribed, candidate species (with statistical 
support from various quantitative species delimitation analyses; 
Fig. 2) exhibit equivalent levels of genetic distinctiveness, at or 
above a level of genetic divergence estimated among formerly 
recognized, names species (Figs. 2, 3). 


Discussion 


Commenting on puddle frogs of northern Borneo, Inger et 
al. (2017) characterized the natural history, ecological 
tendencies, and morphological differences between O. cf. laevis 
and O. baluensis. Describing O. cf. laevis as effectively 
camouflaged against the drab color of the muddy, stagnant water 
in which the species is always found (puddles, small ponds, and 
wallows), the authors described the second species of the 
Borneo dyad (O. baluensis) as matching the variably-colored 
gravel, and red clay soil making up the substate of seeps and 
small flowing streams in forested areas (Table 2). Similarly, our 
survey of genetic, phenotypic, and ecological (microhabitat 
preference) variation in Occidozyga puddle frog populations on 
either side of the Sundaland-Philippine faunal zone interface 
revealed a surprising, widespread pattern of species occurrence 
in dyads, or repeated pairs of distinct evolutionary lineages, 
usually also ecologically constituting somewhat separate “pond” 
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Figure 3. Density plots showing the distribution of uncorrected p- 
distances between reciprocally monophyletic focal clades. 


versus "stream" forms; these involved previously recognized, 
named species, candidate species supported by our statistical 
species delimitation analyses (Fig. 2), and divergent forms with 
phenotypic variation (Table 2), which have diverged from 
congeners on the same scale as divergence found between 
taxonomically uncontroversial named species (Fig. 3). Aside 
from the question of how widespread this pattern may be, and 
whether it prevails on landmasses farther north (the oceanic 
portions of the Philippines) or south (southern Sundaland, and 
other portions of the Indo-Australian Archipelago; Fig. 1), what 
can we infer about evolutionary processes associated with, or 
possibly giving rise to, this newly-elucidated pattern of species 
dyads in geographic distributions of puddle frogs? 

The question articulated above takes a more general form 
if we ask how have recently diverged species pairs (which we 
expect to be phenotypically and ecologically similar, at least at 
the onset of their divergence) empirically been shown to 
coexist? Or, how have seemingly similar species managed to 
coexist in time or space and avoid direct competition? We know 
from numerous comparative phylogenetic analyses of 
ecological niche differentiation, that phenotypically and/or 
ecologically similar land vertebrates appear to coexist due to 
divergence along key axes of organismal variation, often linked 
functionally to fitness via a phenotype-environment correlation 
(Schluter, 2000). In some heavily-cited cases of adaptive 
radiation, primary axes of variation involve evolutionary 
change, in a nearly deterministic fashion, in feeding niche and 
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related specializations (Lamichhaney et al., 2015), specialized 
utilization of structural microhabitats (Mahler et al., 2010), body 
size divergence and existence of optimal size categories (Setiadi 
et al., 2011) or other, critical, form-function relationship features 
(Alfaro et al, 2009; Cooney et al, 2017). In amphibians, 
additional conspicuous potential axes of variation involve larval 
biology (and the entire suite of associated variables, including 
dietary specializations, foraging guilds, life history trait 
variation, developmental timing, etc.; Alcala 1962; Duellman & 
Trueb 1986), reproductive mode, temporal segregation of the 
timing of reproduction, modality of mate-recognition signals 
and acoustic niche partitioning (Duellman and Pyles, 1983; 
Wells, 1977), and the overall fine-scale partitioning of spatial 
differentiation in response to ecological opportunity (Myers and 
Burbrink, 2012; Yoder et al., 2010), which may be tightly 
coupled to fine-scale microhabitat preferences (Blackburn et al., 
2013; Setiadi et al., 2011). Which of these patterns (and inferred 
processes) might explain the species dyad distribution pattern 
revealed here in Philippine puddle frogs? At this point, data are 
quite limited, and these limitations prevent firm conclusions. 
The availability of only a single gene mitochondrial gene locus 
(preventing species tree inferences and characterizations of gene 
flow) and our lack of dense population-level sampling, prevent 
us from confidently distinguishing between geographic structure 
of genetic variation and possible species boundaries (Sukumaran 
and Knowles, 2017); these shortcomings will have to be 
ameliorated with future studies. 

At this point, three primary sources of information hold the 
potential to substantially clarify the questions raised here. First, 
(1) genomic data (Chan et al. 2017; Chan et al. 2020; Hutter et 
al. 2019) collected from throughout the range of our sampling, 
and including key populations on either side of the Sundaland- 
Philippine faunal interface would greatly empower tests of 
species boundaries; second, (2) as the primary anuran mate- 
recognition signal (Wells, 1977), acoustic advertisement calls of 
frogs have the potential to reveal species boundaries, even in 
cases where morphological differentiation has not accompanied 
speciation (Brown et al., 2017, 2015); and third (3) there can be 
no substitute for the kind of natural history information, detailed 
microhabitat data, and behavioral clues that accompany long 
term, comprehensive, survey and re-survey field studies that 
focus on organisms in their natural environment (Brown et al., 
2013b; Sanguila et al., 2016). With their exquisitely derived 
tadpoles (Haas et al., 2014), Occidozyga are known obligate 
carnivores, a fact which opens the possibility of fine-scale larval 
specialization on different prey items in slightly different 
microhabitats (e.g., stagnant pools and puddles, versus flowing 
or cascading, cooler, oxygenated water in forest streams). Thus, 
we strongly encourage students, protected area wardens, wildlife 
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managers, and conservation specialists to join forces in 
collaboration and focus on multilocus genetic datasets (Chan et 
al. 2020), large sample sizes of properly recorded (and 
vouchered) frog vocalizations (Kohler et al., 2017), multivariate 
analyses of adult phenotypic variation (Brown et al., 2017) and, 
most importantly, detailed field-based life history studies of 
Occidozyga microhabitats, larval biology, and adult behavior 
associated with the reproductive effort. Are puddle frog 
populations with apparent preferences for cool, running, 
enclosed forest stream environments truly “specialized” as a 
result of replicated evolutionary shifts in ecological preferences, 
tendencies, or physiological limits—and can we demonstrate 
(with ecological and physiological data and/or diet) that the 
hypothesized, contrasting “pond” forms found in open habitats, 
stagnant pools and puddles, in disturbed areas at low elevations 
have wider water pH tolerances, statistically broader 
temperature range limits, occur in more numerous kinds of 
bodies of water, and/or quantitatively exhibit a more variable 
diet? These and other questions will need to be the subject of 
focal studies, in different areas where puddle frog dyads occur, 
involving populations and species with distinct evolutionary 
histories, and occurring as part of communities of variable 
complexity; real data, collected properly in the field, and 
analyzed with appropriate inferential statistics will be required, 
before vindication or refutation of the pattern hinted at here, 
will be possible. Only then will a taxonomic resolution of this 
complex, archipelago-wide frog group be possible. 

Opening our view through the window of amphibians 
natural history starts with modest, local-scale, observational 
studies requiring patience, persistence, and a desire to 
communicate science (Alcala 1955, 1956, 1957, 1958, 1962; 
Alcala and Brown 1955a,b, 1956, 1982, 1987). Here we have 
shown that rather than representing just a single common, 
widespread, human commensal species, the Philippine 
Occidozyga radiation encapsulates a series of unresolved, 
ongoing, general-interest conceptual questions in ecology and 
evolutionary biology—all of which hold tremendous potential 
for future studies. 
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